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Abstract — Currently, research is being conducted on the optical properties of materials associated 
with the development of solid-state lasers in the 2 micron region. In support of this effort, a mathe- 
matical model describing the energy transfer in a holmium laser sejisitiaed with thulium is developed. 
In this paper, we establish some qualitative properties of the solution of the model, such as non- 
negativity, boundedness, and integrability. A local stability analysis is then performed from which 
conditions for asymptotic stability are obtained. Finally, we report on our numerical analysis of the 
system and how it compares with experimental results. 


INTRODUCTION 

In the early 1980’s, there was a renewal of interest in solid-state lasers, due in large part to 
the development and availability of new host materials. It was ground this time that NASA 
began investigating tunable solid-state lasers as promising candidates for the Earth Observing 
System (Eos). During the past several years, research has been conducted on the sensitized 
holmium laser, which operates in the near infra-red region. It can be used as a source for both 
LIDAR (Light Detection And Ranging) and DIAL (Differential Absorption and Lidar) as well 
as on aircraft to make Doppler lidar measurements of windshear since it operates in the eye-safe 
region. 

In a solid-state laser, a dopant ion substitutes directly into the host lattice. When the lan- 
thanide rare earth ion holmium (Ho) is used as a dopant in the host crystal yttrium aluminum 
garnet (Y 3 AI 5 O 12 ), the holmium ion substitutes into the yttrium sites. This induces a weak 
coupling to the host lattice that results in narrower absorption and emmision features than those 
commonly observed in transition metals. Although the holmium laser has high gain and good 
energy storage properties, ionic interactions among the holmium ions limit the concentration of 
holmium possible in any host; therefore, to increase the optical energy absorbed, a sensitizing 
ion is included. In the laser system under consideration, the lanthanide rare earth ion thulium 
(Tm) is used as a sensitizer for holmium. The interactions between thulium and holmium ions 
increase the efficiency of the Tm-Ho laser but at the cost of introducing more nonlinearities into 
the model. These nonlinearities make both the analysis and the laser dynamics significantly 
different than solid-state laser dynamics studied previously. 

* Author to whom all correspondence should be addressed, 
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The model comprises thulium ions in energy levels 3 H 4 , 3 F<, and 3 H«; holmium ions in energy 
levels 5 I 5 5 It and S I 8 ; and the photon density in the optical cavity. There are two mam types of 
processes’ considered in the laser system; namely, those processes that are inter-ionic-.nvolv.ng 
the transfer of excitation from one ion to another-and those that are not. Inter-iomc processes 
include: cross relaxation, back transfer, and up-conversion. Processes that are not inter-iomc 
include: spontaneous emission, absorption, and stimulated emission. The model does not account 
for spatial dependence of excitation in the crystal rod, but rather utilizes a spacial average o 

the length of the rod. 

In the next section, the system of equations for the electron populations and photon density 
is introduced and discussed. We then establish some qualitative properties of the solution to the 
system After that, a local stability analysis is performed, and finally, the system is subjected 
a numerical treatment and the numerical solution compared with experimental results obtained 

in the laboratory. 


THE MODEL OF THE LASER DYNAMICS 

The energy level diagram for the Tm 3 + ion and the Ho 3 + ion in YAG [1 Figure 2.1, P- 10] is the 
basis for the idealized model, Figure 1, of lasing action For the remainder of tne d,scussm he 
following correspondence is made as a matter of notational convenience: thulium energy levels 
3 H C 3 F 4 and 3 H 4 correspond to energy levels 0, 1, and 2, respectively, and holmium e gy 
levels ft 5 l 7 , and % correspond to energy leve.s O' 1' and 2 , -pecUvely. The numbe^of 
thulium ions, per cm 3 , in energy level * and at time t will be denoted t - 0, 2 ^m.la^ y, 

the number of holmium ions, per cm 3 , in energy level . and at time t will be denoted n,( ), 

i = 0,1,2. 


Pump 
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Figure 1. Idealized model of the Tm-Ho:YAG laser with the follow. ng processes 
represented: (a) stimulated absorption, (b) cross relaxation (c) back transfer, ( ) 
forward transfer, (e) up-conversion, (f) down-conversion, and (g) stimulated emission. 
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The temporal evolution of the dopant electron populations is described by the set of rate 
equations given by 

= W p (t)N 0 (t) - ^ - CN 0 (t)N 2 {t ) 
at r 2 

dNl ^ l = + 2CN 0 (t)N 2 (t) + ClN 0 (t)ni(t) 

dt T21 T\ 

- CiNi(t)no(t) + q[No(t)n 2 (t) - qiN^n^t) 

dN ° ^l = - W p (t)No(t) + - CN Q (t)N 2 (t) 

dt pw v / T2q n 

+ CiN^noit) - ClNo^nxit) + q l N 1 (t)n l (t) - q[N 0 (t)n 2 (t) 

^ M =qiNi(t)ni(t) -?M- q[N Q (t)n 2 (t) 
dt *2 

drnit) = « 2 {t) _ nj^) + ClNl ( t)fl0 ( t) _ ClN 0 (t)ni{t) 
dt r' n t[ 

+ q[N 0 (t)n 2 (t) - qiNi(t)ni(t) - av<l>(t) (n^t) - 

dn °^ l = ^ ^ + ClNoit) ni (t) - CiNi(t)no(t) + avtft) (n^t) - ^■n 0 (t)') - 

dt *20 T i \ 9 0 / 

In this system, C, Cj* , Ci,qi, and q[ represent the probability of cross relaxation, back transfer, 
forward transfer, up-conversion, and down-conversion, respectively. All of these processes are 
depicted in Figure 1. The pumping rate is denoted by W p (t) and represents the number of 
photons, per microsecond, available to excite electrons in energy level 0 to energy level 2. 

The quantities r, , t- i = 1,2 represent the spontaneous emission lifetime of energy level t,t , 
respectively. The quantity l/r,j represents the transition rate (due to spontaneous emission) 
from level i to level j. Similarly, l/r' ; - represents the spontaneous emission transition rate from 
level i' to level j’. These two quantities satisfy the relations 


n 


/=« ■' 


and 


- = Y— 

r f ^ t 1 - • 

T * jzz o »; 


t = 1,2. 


The photon density is denoted by <t>(i), and the rate equation describing the temporal evolution 
of the photon density in the optical cavity is given by 

= ( ni(<) “ ^ no(0 ) _ + sp °^- (2) 

Here, a denotes the transition cross section, v is the velocity of light in the laser crystal, g x is the 

number of manifolds associated with energy level 1', and likewise for 0o* The contribution due 
to spontaneous emission of holmium ions in energy level 3/ is denoted sp Cl and the fluorescent 
lifetime of the material is denoted ryj. The lifetime of a photon in the optical cavity is denoted 
by r c and can be expressed as [2] 

- 2 ^/ c 
Tc ” — In RiRn ’ 

where C c is the length of the optical cavity, c is the velocity of light in vacuo , and R x and R 7 
denote the reflectivity of the two mirrors. The electron populations of thulium and holmium are 
constrained by the relations 

Nt = No(t) + N\(t) + A^(0 and nr = «o(0 + n i(0 + n 2(0» 0) 

where N T is the concentration of thulium ions, per cm 3 , and n T is the concentration of holmium 

ions, per cm 3 . 


KM 17.6-f 
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These “population constraints” can be rewritten as 

N 0 (t) = N T - Ni(t) - N 2 (t) and n o (0 = "T - «i(0 “ "2(0- ( 4 ) 

Substituting these expressions for Wo(0 and n 0 (t) into system (1) reduces the system from six 
to four equations. Then, by including the rate equation for the photon density (2), we obtain 


- w p (t)(N T - Ni(t) - N 2 (t)) - - CN 2 (t)(N T - Ni{t) - N 2 (t)) 

dt 

MW) = WA _ *W) + 2 CN 2 (t){N T - Ni(t) - N 2 (t)) 
dt 7*21 n 

+ <?i*"i(0(Wr - Wi(0 - Wz(O) - CiNi(t){n T - m (t) - n 2 (t)) 

+ q'lWWr - Nx(t) - N 2 (t)) - 

= qi N\(t)ni(t) - ^ - q[n 2 (t)(N T - Ni(t) ~ N 2 (t)) 
dt T 2 ' * 

MWt _ _ HM + CiNi(t)(nT - ni(t) - n 2 (0) 

dt ^2i 'i 

- Clni(t)(N T - Ni(t) - N 2 (t)) + q[n 2 (t)(N T - Ni{t) - N 2 {t)) 

- tMt)ni(t) - <rv<!>{t) ^m(<) - - "i(0 “ "2(0)) 

tm = (».«) - £<■.«)) - ^ + **^r- 

The coupled system (5) is the temporal model of the laser system expressed in terms of the 
physical variables. We now normalize system (5) to put it into the form used for the analysis 
and calculations in the subsequent sections. Hence, let 


x(0 = 


Ml 

Nt * 


y(0 = 


Ni(t) 
Nt 5 


z(0 = 


"i(0 

i 

rtT 


w(t) = 


" 2(0 

I 

7lT 


P( 0 = 


Hi] 

^norm 


( 6 ) 


be the normalized variables. Rewriting the rate equations (5) in terms of the normalized vari- 
ables (6) yields 


where 


% = W„(0(1 - x - y) - — - D lX (l - x - y) 
dt t 2 

+ 2L»ix(l - x - y) + D 3 z(l - x - y) 

dt X2i n 

- D 2 y{\ - z - w) + D 6 iy(l - x - y) - D 7 yz 

— D *v z ~ ~j ~ ^ 9UI (i ~ x ~y) 

dt *2 

± = ^--± + D 4 y(l -z-w)- D s z(l - x - y) 
dt Tj X tJ 

+ D 9 u>( 1 — x — y) - D 6 yz - P\[~f z + (1 - t)(1 — W )\P 
^ = jft[7* + (l-7)(l-«')]-^} P + 03* 


Di = CN t 

D 2 = CitlT 
D 3 = C[nx 



Da = C\Nt 

d 5 = c;nt 

Dq = {jnr 


£> 7 = qiTlr 
Ds = qiNr 
Ds = q[Nr 


Pi — <rv^ n0 rm 
p 2 = a-imr 
. spo n T 

P3 = 


7// 


(7) 
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QUALITATIVE ANALYSIS 

In this section, the full set of equations (7) is modified by excluding two physical processes— 
up-conversion and back transfer. This yields a simplified system for which certain qualitative 
properties of the solution are established. Throughout this discussion, the parameters are as- 
sumed to satisfy the following (physically realistic) relationships: 


r 2 1 > T2 > 0 
Tj > r 2 o > 0 

Tj > t c > 0 

I @2, 03 > 0 

Di,Di,Di > 0 

K 7 


Sp 0 << 1 . 

Excluding the above mentioned processes yields the following simplified system of rate equations 


= W,(t ) (1 - x - y)+ - ~ - L>ix(l - * ~ 1/)+ 


where 


^ = — - 2- + 2D 1 x(l - x - y) + - D 2 y( 1 - z)+ 

dt T 2 1 Ti 

^ = --7 + D 4 y ( I - z)+ + /?i(7(l “ *)+ ~ 

dt T[ 

~ = j&[l “ 7(1 “ 2 )+] ~ P + 



1 — x — y, if 1 — a: — y > 0, 
0, otherwise, 

1 - z, if 1 - z > 0, 

0, otherwise. 


( 8 ) 


(9) 


It follows directly from [3] that if (i) W p (t) is continuous, and (ii) the solution vector, y - 
[x,y,z,P] T , is bounded, then there exists a unique, continuous solution to the system (8). The 
following four theorems may be viewed as a validation of the model since they show that the 
solution to system (8) possesses the same kind of qualitative properties that one would expect 
the laser system to possess. Selected portions of the proofs are given below. The portions that 
are omitted can be found in [1, pp. 48—68], The first theorem establishes the nonnegativity of 
the solution vector y, as well as the ground states of thulium and holmium. 


Theorem 1. (Nonnegativity). UW p (t) > 0 for t > 0 and x(0) > 0, y(0) > 0, *(0) > 0, and 
P(0) > 0, then 

(i) If 1 - x(0) - y(0) > 0, then 1 - x(t) - y(t) > 0 for all t > 0; 

(ii) If 1 — x(0) - y(0) < 0, then there exists a T such that 1 — x(f) - y(<) > 0 for all t > T; 

(iii) (a) If 1 - x(0) - y(0) > 0, then x(t) > 0 and y(<) > 0 for all t > 0 
(b) If 1 - x(0) - y(0) < 0, then x(t) > 0 and y(t) >0for allt> 0; 

(iv) If I - z(0) > 0, then 1 - z(t) > 0 for all t > 0; 

(v) If 1 - z(0) < 0, then there exists a T* such that 1 - z(t) > 0 for all t > T* ; 

(vi) z(t) > 0 and P(t) >0 for all t > 0. 


Proof of Parts (i), (ii), and (iii) of Theorem I 

Proof (i). Since z(0) + y(0) < 1, then by continuity ofx(t) + y(f), we know that i(<) + y(t) < 1 
for t near zero. Suppose that at t = T, x(T) + y(T) = 1 and x(t) + y(t) < 1 for 0 < t < T, then 


(10) 
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Now since l/r 21 - l/r 2 < 0 and -1/n - D 2 ( 1 - z) + < 0, then if x(T) > 0 and y(T) > 0 and not 
both are zero, then from (10), £(x + y) < 0 at t = T which contradicts x(<) + y{t) < 1 for t < T 
and therefore x(<) + y{t) < 1 for all t. Now since x(T) + y(T) = 1, not both x and y are zero at 
t = T\ hence, it remains to be shown that x(T) > 0 and y(T) >0. 

Suppose y(T) < 0, then since y(0) > 0 and y is continuous, there exists a T y <T such that 

y(T v ) = 0 and y(<) > 0 for 0 < t < T y . Then at t = T y 

— + 2L>i(l -x-y)+ . 

T 21 J 

The quantity in brackets above is positive, so if x(T y ) > 0, then y(T y ) > 0 which means y is 
increasing at i — T y . But this contradicts y(t) > 0 for 0 < t < T y , and this contradiction implies 
y(<) > 0 for all t. Now suppose x(T y ) < 0, then since x(0) > 0 and x is continuous, there exists 
a T x < T y such that x(T*) = 0 and x(t) > 0 for 0 < t < T x . Then at t = T x 

^ = W P (T x )(l - x - y) + > 0. 

Hence, x(T x ) > 0, which means x is increasing at i = T x , which contradicts x(t) > 0 for 
0 < t < T x . Therefore, x(T y ) > 0 and x(t) > 0 for 0 < t < T y . Now suppose x(T v ) = 0, then 
y(T y ) = 0 and x(T y ) = W p {T y ) > 0. Additionally, 

+ 2£) 1 x-j-(l - x - y) + + 2£>ix(l - x - y) + 
dt 2 Toi T\ dt 

- Diy J t (l - z)+ - D 2 y(l - z) + , 

from which we obtain 

g (71)= (j r+2A ) i(J; ) 

> o. 

So y{T y ) = 0, y{T y ) = 0, and y(T y ) > 0. Hence, y is concave up at t = T y , which implies that y 
is increasing to the right of / = T v , and so 

y(T) > 0. (11) 

Now suppose x(T) < 0, then since x(0) > 0 and x is continuous, we can choose T) < T such that 

x(T\) = 0 and x(<) > 0 for 0 < t < Ti, then at t = Ti 

^ = W p {T l )( l-x — y)+ 

>0, 

so x is increasing at t = Ti which contradicts x(t) > 0 for 0 < / < Tj ; hence, the above supposition 
is false which implies x(T) > 0. This, together with (11), shows that at t = T 

7t {x + y)<0 - 

From prior comments, this implies that x(t) + y(t) < 1 for ail t . I 

Proof (ii). The proof of (i) shows that if x(Q) + y( 0) = 1, then x(t) + y(t) < 1 for all f > 0. 

Suppose x(0) + y(0) > 1, then for t near zero 

7i (x * s) = H + * ' + ('n ■ m ' l)+ ) 51 

<0, 
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since for t near zero -l/r 2 + l/r 21 < 0, -1/n - D 2 { 1 - z) + < 0, x(t) > 0, y(t) > 0, and not 
both are zero at t = 0. So x(t) + y(t) is decreasing for t near zero. Now assume there exists a £ 

such that x(t) + y(<) > ^ > 1 for all t, then 

= -i- + — ^ x — — (x + y) — D 2 (l — x)+y 

V T2 T21 n/ r i 

<--C 

T-l 

, /_ , 1 . i / T , <n —Di(l — x)+ < 0, and an argument similar to that given in (i) 

“"e iduLto show that a' > 0, y > 0. This inequality holds for all I, so upon integral^ 
from 0 to t, we obtain 


x(0 + y(t) - (x(0) + y(0)) < 


or equivalently, 

From inequality (12), and since n > 0, it follows that 


x(*) + y(0 < & + x (0) + 


( 12 ) 


lim (x(<) + y(0) = -°°i 
<—►00 


which contradicts *(!) + „(<) > « > 1 for all t. Hence our assumptron “P 1 " 

there exists a T such that x(T) + »(T) = 1 and *(t) + !>(<) <‘ f °> ‘ "ear T. The analysis from (r) 
now applies to show that x(<) + y(<) < 1 for all t > T as required. ■ 

PROOF (iii(a)) . Let x(0) + y(0) < 1, then from (i), we know x(t) + y(t) < 1 for all t. Also by 
hypothesis, ^ 7 ( 0 } > 0. Therefore, since x is continuous, x(t) > 0 for t near_zero. Assume at some 
point, say t = T x , that z(Ti) = 0 and x(t) > 0 for 0 < t < Ty Then at t - T u 


^ = WV(T0(1 - y(Ti))+ 
at 

' > 0 , 


which means x is increasing. This contradicts x(t) > 0 for 0 < t < Tj. Hence, the above 
assumption is false, which implies x(t) > 0 for all t > 0. 

Similarly, assume that at some point, say t = T 2 , y(T 2 ) = 0 and y(t) > 0 for 0 t _ 2- 
at t = T 2 , we have 


positive 

^ = f— + 2Di(l — x)+l x 
dt [t 2 i J 

> 0 , 


since x > 0 and l/(r 21 ) + 20^1 - x) + > 0. This implies that 
contradicts y(t) > 0 for 0 < t < T 2 . Hence, the assumption was 


y is increasing at t = T 2 which 
fallacious which shows y(0 > 0 


for all t > 0. 

Proof (iii(b)). Let x(0)+y(0) > 1, then from (ii) we know there 
1 for all t > T, and x(T) + y(T) = 1. Suppose that at T, < 
0 < t < 7\ . Then at t = T x 

— — ri r 


exists a T such that x(t)+y(0 < 
T, x(Ti) = 0 and x(t) > 0 for 
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Also, at t = T\ we have 


^(T t ) = -^ p (T 1 )y(T 1 ) 


= - ^>2sK^i)(i - *)+ 

= W p (Ti)y(Ti) Ji + D 2 (l-z) + ] 


positive 


>o. 


Hence, x is concave up at t = T\ and consequently increases to the right of it. From this, it 
follows that x(<) > 0 for 0 < t < T\ furthermore, for t > T the analysis from (iii(a)) applies to 
show x(t) > 0 for all t > T. 

Again, assume that at T 2 < T t y(T 2 ) = 0 and y(t) > 0 for 0 < t < T 2} then at t = T 2 


dy 

dt 


±-x(T 2 ) 


> 0 . 


So y is increasing at t = T 2 which contradicts y(t) > 0 for 0 < t < T 2 ; hence, the assumption was 
fake which shows that y(t) > 0 for 0 < t < T. Ako, for t > T, the analysis from (iii(a)) pertains 
to show y(t) > 0 for t > T. Taken together, we have y(t) > 0 for all i > 0. I 

The proofs for parts (iv), (v), and (vi) proceed along similar lines. 

Theorem 1 states that the normalized ground levels of thulium and holmium remain nonneg- 
ative on I = [0,oo), or, in other words, that the constraints (9) are self-enforcing. Consequently, 
we will drop the “+" subscripts in system (8) for the remainder of the discussion. 

Theorem 2. ^Integr ability). Let W p (t) be positive, continuous, and integrable on I. Then 
ifx( 0) > 0, j/(0) > 0, z(0) > 0, and P(0) > 0, then x(t), y(t), z(t), and P(t) are all integrable 
on I. 

PROOF. Fix T such that x(t) -I- y(t) < 1 for all t > T, then for t > T we have 

$ = -- - D lX ( 1 -x-y)+ W p (t) - W P (t)x - W P (t)y 
at t 2 

= (“ - ^p(o) 1 + w tW - A*(i - x - y) - w)>(0y 

<-(y 2 + W p (t)jx + W p (t), 

since W p (t) > 0, y > 0, (1 - x - y) > 0, and x > 0. Rearranging, we get 

Now multiply both sides by E(t) = t j T (( l l T ?)+ w r(0) d ( 0 b^aj n 

£(0§ + (j- + W,(t)) E(t)x < E(t)W p (t), 


which implies 


^{£(t)*(<)} < E{t)W p {t). 


(13) 
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Upon integrating inequality (13) from T — * t and noting that E(T) — 1, we obtain 

E(t)x(t) - x(T) < [ W p (u)E(u) du, 

JT 

which holds for t > T. It now follows that 

0 <*(<)< [ WpMEME-'Wdu + xiDE-'it). 

J T 

Since 


(14) 


can express (14) as 

0 <*(<)< ! '- W ' MW ' 


Jt 

5 , integrating from T — ♦ s yields 




, UltC£i nuiu j. r « j 

f x(t)dl < jf |jf’ iv,(») e - i*} * + 

iince tVp is continuous on I, the order of integration may be interchanged to obtain the 
ing equivalent expression 

f x(t)dt < /; { jf' *} * + «T) jf 

the relation 


dt. 

(15) 


e 


- f* JV.fCWC 


it follows from (15) that 

J ’ x(t)dt < j r ^p(u) {jf dt} du + x(T) {jf e- ( ‘- T)/ra *} 

= J’ t 2 W p (u){\ -e-(*- u)/Ta } du + x(T)r 2 {l-e- (, - T)/T3 }. 

Since r 2 > 0 and Wp(f) is integrable on X, it now follows that x(t) is integrable on X. 

To show that y(t) is integrable on J, start with 

— (x + y) = ( — — + — ^ x - — — £> 2 y( 1 — z) + D\x(l - x — y) + W p — W p x - W p y 
dt' ' V r 2 T 2i/ n 

and proceed along the same general lines as above to establish that (x + y) is integrable on X. 
Then invoke the nonnegativity result for x and y to conclude that y(t) is integrable on 1. 

Theorem 3. (Boundedness). If W p {t) > 0 on X, then x(t), y(t), z(t), and P(t) are a II bounded 
onl. 

Theorem 4. Let W p , x, y, z, and P satisfy the hypotheses of Theorem 1. Furthermore, let 
W P (t) —* 0 as t — ► oo. Then, the solution vector y — > 0 as t — ► oo. 
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STABILITY ANALYSIS 

In this section, the equilibrium solutions are obtained, and a local stability analysis is per- 
formed. Before doing this, however, we make the simplifying assumption that the spontaneous 
emission term, sp 0 , has a negligible affect on the asymptotic behaviour of the system and hence 
may be dropped. Furthermore, we will consider the case of continuous wave (CW) pumping. 
These two considerations are tantamount to taking ft = 0 and Wp(t) = W P in system (8), where 
W p is a positive constant. Doing this yields the following system of rate equations: 

^ = W p ( 1 — x — y) ~ Dix(l - x-y ) 

— = — — — + 2DiJc(l — x — y) — D 2 y(l — z) 

dt t 2 i n ( 16 ) 

^ = -I T + D 4 y(l-z) + P 1 [r(l-z)-l]P 

at T* 

f = {ftP -70-)] -£}« 

By setting the right hand side of equations (16) equal to zero, two equilibrium points are obtained. 
Equilibrium point 1 is denoted I, which has coordinates . (x 1 , y ,z ,0), and equilibrium point 2 
is denoted II, which has coordinates (z 2 , y 2 , z 2 , P 2 ), where the coordinates are given in terms of 

the physical parameters of the system. 

Having located the equilibrium points, we now determine the stability properties of I and II. 
For the stability of I, introduce the new hat variables given by x = x-x l i y — y-y l ) z^z-z , 
and P = P - 0. Linearizing system (16) about I yields 

i = + ( 17 ) 

where A represents the coefficient matrix of the linearized system, X = [x } y,z t P] T } and the g 
vector contains all the nonlinearities of the system. 

Invoking the Routh-Hurwitz Theorem [1, pp. 71-73], the following necessary and sufficient 
condition for I to be asymptotically stable is obtained: 

ft[l-7(l-^)]-7<°- 

*C 

The stability analysis for II parallels that of I and yields one condition for II to be asymptotically 
stable. The algebra associated with this condition, however, is tremendously complicated, and 
hence, we resort to a parametric investigation of the stability of the equilibrium points. 

The numerics indicate there is an interplay of stability between I and II. To see this, we 
consider W p as a parameter and choose it to be “small ” Then using a computer program, we 
calculate the equilibrium values of I and II and use the Routh-Hurwitz criteria to determine the 
stability of each point. By gradually increasing W p and repeating the process, we see there is an 
interchange of stability which goes like this: for W p small, the P coordinate of II is negative, II 
is unstable, and I is asymptotically stable (i.e., no lasing). As W p is increased, the P coordinate 
of II eventually reaches zero, and thus, the two equilibrium points coalesce (this is the threshold 
of lasing). Finally, as W p is increased still further, the P coordinate of II becomes positive, II 
becomes asymptotically stable, and I becomes unstable (i.e., lasing occurs). 

This is illustrated in Figure 2 where P 2 is graphed as a function of W p , and the stability of I 
and II is noted. The value of W p at which the interchange of stability occurs is a bifurcation 
point and is denoted W*. From numerical calculations, we find W* = 0.342 x 10 . 

It is also of physical interest to see how the system responds as the cavity lifetime of a photon, r Cl 
is varied. Choosing r c small and following the same procedure as above yields the same interplay 
of stability as before. Figure 3 gives P 2 as a function of t c , once again noting the stability of 
each point. From this, we see that r c has a bifurcation point which will be denoted r* and is 
approximated by r* =4.81 x 10 -5 . 
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Figure 2. P coordinate of II as a function of the pumping term. 



Figure 3. P coordinate of II as a function of the cavity lifetime r c . 


NUMERICAL ANALYSIS 

In this section, the numerical solution to system (8) is obtained. This is done by employing 
either the subroutine LSODA (when using a VAX 11/750) or the subroutine DDRIV2 (when 
using an IBM PC). In either case, all computations were performed in double precision arith-. 
metic. LSODA was developed in 1987 at Lawrence Livermore National Laboratory in Livermore, 
California, by L. R. Petzold and A. C. Hindmarsh. DDRIV2 was developed in 1979 and re- 
vised in 1987 by D. K. Kahaner, National Bureau of Standards, and C. D. Sutherland of Los 
Alamos National Laboratory. Both subroutines were created for the numerical integration of stiff 
and nonstiff systems of first-order ordinary differential equations. When the system is stiff, the 
subroutines use a backward difference formula (BDF) to perform the numerical integration, and 
when the system is not stiff, they use a higher-order Adams method. Both LSODA and DDRIV2 
have the capability of automatically switching from one method to the other as the system passes 
from stiff to nonstiff regions. The BDF method was chosen by both d.e. solvers throughout the 
interval of integration which indicates that system (8) is stiff. Roughly speaking, a system is 




78 


T.G. WANGLER ct at. 


stiff when the characteristic equation of the associated linear system has a root with a “large” 
negative real part. A more detailed treatment of stiffness can be found in [4]. The numerical 
values of the system parameters are given in Table 1 and were used in the computation of the 
numerical solution. Employing the program in [1, Appendix B], the numerical solution of (8) was 
obtained. A typical example of the numerical solution is given in Figure 4 where the upper lasing 
level and the photon density are plotted. Since W p is constant, a comparison of the Tong time” 
behaviour of the numerical solutions with that of the equilibrium solutions can be performed. 
Using the parameter values in Table 1 and the Routh-Hurwitz criteria for asymptotic stability, it 
is found that II is asymptotically stable and I is unstable. The coordinates of II, as predicted by 
the stability analysis, are given in the left column below, and the numerical solutions at t = 200/is 
are given in the right column below. 


X 2 = 0.1304363 

x = 0.1304372 

y 2 = 0.7637408 

y =0.7637406 

z 2 = 0.5238095 

z =0.5238095 

P 2 = 0.0759492 

P = 0.0759491 


From this, it is evident that the d.e. solver captures the Tong time” behaviour of the solutions 
extremely well. 

Table 1. System parameters for Tm-Ho:YAG Laser. 


N t = lx 10 21 /cm 3 

T2i = 900 /is 

Tir = 1 X 10 2 °/cm 3 

ti = 11,000/iS 

t = 0.3 cm 

Tj = 8, 500 /is 

<£norm = 1 X 10 /cm 

$ 

II 

»-« 

t c - 17.5 cm 

9l = 1 

R x - 1.0 

7 = 2 

R 2 = 0.95 

a = 7 x 10“ 21 cm 2 

C = — — — cm 3 /(is 

o 

6 

II 

51 

40 x 10 21 

1 3 , 

It 

O 

o 

Cx = — cm J / /is 

475 X 10 20 

v = 30, 000 cm/ /zs 

o 
*“* • 

II 

o 

o 

T C = 1 X 10 ” 3 /iS 

T2 = 450 /IS 

W p = 6 X 10“ 3 /cm 3 • fis 

T 20 — 900 /is 

spo — lx 10~ a 


Various pumping schemes were considered in the model. Figure 5 gives the numerical solution 
for P(t) when the pumping term is taken to be 

W p {t) = a 2 <e-“ V , 

where a = .01133. From this figure, we see that the photon density has irregular oscillations and 
decays rather rapidly, becoming negligible by about t = 250/rs. 

We now consider two modifications that are made to upgrade the model. The first is to include 
the terms that account for the back transfer of energy and that were excluded at the beginning 
of the section on qualitative analysis. The second is to replace 02 and 0 3 in the rate equation for 
the photon density with 

02 = (j-j 02 and 03 = 03. 

The rationale for the second change is as follows. A mathematical model that accurately describes 
the dynamics of a laser system must account for both spatial and temporal variations in the 
dependent variables. Allowing for both types of variation yields a system of nonlinear partial 
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CONSTANT PUMP (WP = .006) 


Figure 4 . Normalized photon density and upper lasing level with a constant pump. 
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Figure 5. Normalized photon density with a decaying exponential pump. 

differential equations (p.d.e.*). Due to the intractability of the equations and the amount of effort 
and computer time involved in solving the system numerically, various averaging techniques are 
used to eliminate the spatial dependence and thus make the system more tractable. L. F. Roberts 
ei ai [5] gives three such averaging schemes along with their resultant o.d.e. models. The con- 
clusion of this work is that only the temporal model obtained by taking a spatial average over 
the optical length of the cavity, l Cy yields numerical results that agree qualitatively with those 
predicted by the spatial and temporal model. This averaging scheme can be incorporated into 
the model by replacing /?2 and /? 3 by the quantities given above. Doing this and including the 
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back transfer terms yields the system 

^ = W„(0(1 -x-y) - — - D x x{\ - x-y ) 
at T2 

= — - 2- + 2D 1 x(l -x-y)- D 2 y{ 1 - z) + D 3 z{ 1 - x - y) . 
dt t 2 j n (i8) 

— = —L + £> 4 y( 1 -z)- £) s ^(l - x - y) + - *) - 1]** 

dt tJ 

The numerical solution to the system (18) was obtained giving special attention to the qualitative 
behaviour of the photon density for various values of Cl and t/t c - Recall that C* represents the 
probability of back transfer occurring. The pumping term was taken to be of the form 

W p (t) = a 2 ^- 0 ’* 3 . 

In Figure 6, the normalized photon density is graphed for a set of parameter values that are 
“close” to the values used in experimentation. Comparing Figure 6 with Figure 5, we see a 
salient disparity in the qualitative behaviour of the photon density. The previously observed 
erratic spiking has been replaced by regular and temperate oscillations. Figure 7 is a picture of 
the energy output of the laser system as displayed on the screen of an oscilloscope. By comparing 
Figure 6 with Figure 7, we see that the behaviour of the photon density, as predicted by the 
temporal model (18), is in remarkable agreement with the behaviour observed in the laboratory 
(both having about 7 oscillations in 25 /is). 



Figure 6. Photon density when using an alternate averaging technique. 

Finally we further upgrade the model by including the terms associated with the up-conversion 
of energy’ from the 5 I 7 energy level to the 5 I 5 energy level of holmium. This now brings us full- 
circle since we are back to the full set of equations (7) with the exception that ft and ft are 
replaced by (3' 2 and 0' 3l respectively, as defined above. The numerical solution to the system with 
up-conversion was computed taking the probability of up-conversion and down-conversion as 

qi = 5.0 x 10~ 22 and q[ = 1.2 x 10 -23 , 

respectively. Figure 8 shows the interplay between the photon density and the upper lasing level 
of holmium, and Figure 9 gives a phase portrait of the same. 
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Figure 7. The energy output of the Tm-Ho:YAG laser as displayed on the screen of 
an oscilloscope. 
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Figure 8. Upper lasing level and photon density with up-con version included in the 
model. 


SUMMARY AND CONCLUSIONS 

We have developed a temporal model of the dynamics of an optically pumped co-doped six-level 
solid state laser. The model was developed to study the inter-ionic transfer of energy between 
thulium and holmium and how this affects the performance of the laser system. In this paper, 
we reported on the quantitative and qualitative behaviour of the solutions. 

The validity of the simplified model was established by showing that — under appropriate con- 
ditions on the pumping term — solutions to the system exist and are bounded. Further, we showed 
that if the initial conditions are physically realistic and the pumping term is well-enough behaved, 
then the solutions are both nonnegative find integrable (integrability is desirable since calculating 
the efficiency of the laser system involves integrating over the photon density). 

Two equilibrium solutions were obtained and a local stability analysis performed. The interplay 
of stability between the two equilibrium points was then demonstrated parametrically. 

The system was solved numerically and found to be stiff. The d.e. solvers used performed 
exceptionally well and gave numerical solutions that agreed with the stability analysis and the 
experimental results. By gradually upgrading the simplified model, it was found that both back 
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Figure 9. Phase portrait of the upper lasing level and the photon density. 

transfer and up-conversion affect the time of lasing as well as the magnitude and frequency of 
the laser output, but it is an alternate averaging scheme that has the most dramatic affect on 
the qualitative behaviour of the output pulse. 
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